#include "AnaDef.h"

class AnaTPCdEdxHist
{
public:
  static Double_t Ymin(){return 0;}
  static Double_t Ymax(){return 200;}
  static Double_t Ny(){return 400;}
  static Int_t Nx(){return N1x()+N2x();}
  static Double_t Xmin(){return -1; }
  static Double_t Xmax(){return 2; }
  static void SetPPtAxis(TAxis *ax);

private:
  //< Xmid higher bin density
  static Double_t Xmid(){return 0.6;}
  static Int_t HigherDensity(){return 10;}
  static Int_t ScaleUpNXbins(){return 1;}//<--- test
  static Int_t N1x(){return 32 * HigherDensity() * ScaleUpNXbins();}//8*4 * 10
  static Int_t N2x(){return 28 * ScaleUpNXbins();}//4*7
};
  
void AnaTPCdEdxHist::SetPPtAxis(TAxis *ax)
{
  const Double_t dx1 = (Xmid()-Xmin())/N1x();
  const Double_t dx2 = (Xmax()-Xmid())/N2x();

  Double_t bins[Nx()+1];
  bins[0] = Xmin();
  //[Xmin, Xmid] N1x bins
  for(Int_t ibin=1; ibin<=N1x(); ibin++){
    bins[ibin] = bins[ibin-1]+dx1;
  }
  //[Xmid, Xmax] N2x bins
  for(Int_t ibin=1; ibin<=N2x(); ibin++){
    bins[N1x()+ibin] = bins[N1x()+ibin-1]+dx2;
  }

  ax->Set(Nx(), bins);

  /*
  printf("\n");
  for(Int_t ibin=0; ibin<=ax->GetNbins(); ibin++){
    const Double_t xh = ax->GetBinUpEdge(ibin);
    printf("AnaTPCdEdxHist::SetPPtAxis %d: %e = log10 %e\n", ibin, xh, TMath::Power(10, xh));
  }
  printf("\n");
  */
}

int main()
{
  //================================================================
  //                           initialization
  //================================================================

  const Int_t nvar = 5;
  const TString atitle[nvar]={"eta","Pt_jet", "TPCsig", "Log10P_track", "Log10Pt_track"};
  const Int_t nbins[nvar]   ={4,    4,   AnaTPCdEdxHist::Ny(),   AnaTPCdEdxHist::Nx(),                  10000};
  const Double_t xmins[nvar]={0,   20, AnaTPCdEdxHist::Ymin(), AnaTPCdEdxHist::Xmin(), AnaTPCdEdxHist::Xmin()};
  const Double_t xmaxs[nvar]={0.9, 80, AnaTPCdEdxHist::Ymax(), AnaTPCdEdxHist::Xmax(), AnaTPCdEdxHist::Xmax()};

  THnSparseD * hh = new THnSparseD("anatpcdedxhist","", nvar, nbins, xmins, xmaxs);
  for(Int_t ia=0; ia<nvar; ia++){
    hh->GetAxis(ia)->SetTitle(atitle[ia]);
  }

  TAxis * ax = 0x0;
  Int_t nb = 0;

  //0: eta bin
  ax = hh->GetAxis(0);
  nb = ax->GetNbins();
  const Double_t etas[]={0, 0.2, 0.4, 0.6, 0.9};
  ax->Set(nb, etas);

  //1: Pt_jet bin
  ax = hh->GetAxis(1);
  nb = ax->GetNbins();
  const Double_t ptjets[]={20, 30, 40, 60, 80};
  ax->Set(nb, ptjets);

  //2: TPCsig bin equal width, no need to set

  //3: Log10P_track bin
  AnaTPCdEdxHist::SetPPtAxis(hh->GetAxis(3));
  
  //4: Log10Pt_track bin equal width, no need to set

  //================================================================
  //                           test
  //================================================================
  TRandom3 ran(1);
  for(Int_t ii=0; ii<10000; ii++){
    const Double_t eta = -0.9 + 1.8*ran.Rndm();
    const Double_t jetpt = 20 + 60*ran.Rndm();
    const Double_t tpc = 200*ran.Rndm();
    const Double_t logp = -1 + 3*ran.Rndm();
    const Double_t vars[]={eta, jetpt, tpc, logp, logp};

    hh->Fill(vars);
  }

  TFile *fout=new TFile("outtest.root","recreate");
  hh->Write();
  fout->Save();
  fout->Close();

  delete fout;

  return 0;
}
